Hierarchical Time Series Models

library(data.table)
library(tidyverse)
library(lubridate)
library(xts)
library(leaflet)
library(ggfortify)
library(tidyquant)
library(forecast)
library(tseries)
library(TSA)
library(hts)
library(fpp2)
library(tsbox)
rm(list = ls())
gc()

##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 2289029 122.3    4170450 222.8         NA  4170450 222.8
## Vcells 3858122  29.5    8388608  64.0      16384  7003851  53.5

Import Data

datapath = "/Volumes/SSD/TS-Project/data"
train = fread(file=paste(datapath,"subset/train_city_type_item.csv",sep="/"))

Transform to Wide form

# Long to Wide
train_wide = spread(train, hts_label, tot_sales)
train_wide[is.na(train_wide)] = 0
train_wide = train_wide %>%
                  select(-V1, -city, -type, -family, -price, -n_promotion, -perishable, -national_flag,
                         -holiday_flag, -pay.day_flag) %>%
                  group_by(date) %>%
                  summarise_all(funs(sum))
# Covert to TS object
train.ts = ts(train_wide[2:ncol(train_wide)], frequency=365.25, 
     start = c(year("2013-01-01"), 1))

Run HTS

train.hts = hts(train.ts, characters = c(2, 1, 5))
# summary(train.hts)
# smatrix(train.hts)

Plot HTS

# Graph by City
train.hts %>% aggts(levels=0:1) %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City")

# Graph by Type - City Guayaquil
aggts(train.hts, levels=1:2)[,c(1,3:7)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City & Type - Guayaquil")

# Graph by Type - City Quito
aggts(train.hts, levels=1:2)[,c(2,8:11)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City & Type - Quito")

# Graph by Item Family - Type A, City Guayaquil
aggts(train.hts, levels=2:3)[,c(1,10:14)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Guayaquil Type A")

# Graph by Item Family - Type B, City Guayaquil
aggts(train.hts, levels=2:3)[,c(2,15:19)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Guayaquil Type B")

# Graph by Item Family - Type C, City Guayaquil
aggts(train.hts, levels=2:3)[,c(3,20:24)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Guayaquil Type C")

# Graph by Item Family - Type D, City Guayaquil
aggts(train.hts, levels=2:3)[,c(4,25:29)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Guayaquil Type D")

# Graph by Item Family - Type E, City Guayaquil
aggts(train.hts, levels=2:3)[,c(5,30:34)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Guayaquil Type E")

# Graph by Item Family - Type A, City Quito
aggts(train.hts, levels=2:3)[,c(6,35:39)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Quito Type A")

# Graph by Item Family - Type B, City Quito
aggts(train.hts, levels=2:3)[,c(7,40:44)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Quito Type B")

# Graph by Item Family - Type C, City Quito
aggts(train.hts, levels=2:3)[,c(8,45:49)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Quito Type C")

# Graph by Item Family - Type D, City Quito
aggts(train.hts, levels=2:3)[,c(9,50:54)] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Grocery Sales by City, Type & Item Family - Quito Type D")

Train - Val Split

data.hts = window(train.hts, start = 2013, end = 2017.541) # 1659 dates
summary(data.hts)

## Hierarchical Time Series 
## 4 Levels 
## Number of nodes at each level: 1 2 9 45 
## Total number of series: 57 
## Number of observations per series: 1659 
## Top level series: 
## Time Series:
## Start = 2013 
## End = 2017.53935660506 
## Frequency = 365.25 
##    [1] 241577.0 171324.0 164473.0 235336.0 265571.0 160292.0 153524.0
##    [8] 145601.0 121120.0 133692.0 189548.0 233075.0 138496.0 137867.0
##   [15] 158217.0 127061.0 139475.0 208187.0 232920.0 146096.0 145599.0
##   [22] 135308.0 112211.0 133638.0 196931.0 213208.0 130887.0 122364.0
##   [29] 136856.0 123932.0 172838.0 246361.0 242835.0 154577.0 147509.0
##   [36] 150533.0 130835.0 154613.0 177656.0 133130.0 129093.0 153207.0
##   [43] 143646.0 147843.0 147110.0 240149.0 197922.0 168330.0 135163.0
##   [50] 151179.0 121412.0 129113.0 195541.0 226900.0 142030.0 123816.0
##   [57] 135900.0 127229.0 173767.0 263435.0 285592.0 179279.0 158242.0
##   [64] 161989.0 135163.0 149503.0 217245.0 232390.0 147352.0 125691.0
##   [71] 140955.0 112799.0 146067.0 235891.0 265151.0 153670.0 136522.0
##   [78] 157052.0 124680.0 139116.0 220968.0 227613.0 148319.0 126721.0
##   [85] 153239.0 160212.0 167735.0 214670.0 237042.0 201574.0 177582.0
##   [92] 179609.0 141698.0 153072.0 248134.0 264936.0 161229.0 138056.0
##   [99] 149411.0 141702.0 143868.0 219613.0 243008.0 160746.0 160701.0
##  [106] 159686.0 126704.0 142720.0 222631.0 249287.0 154733.0 139622.0
##  [113] 150895.0 125275.0 139525.0 213272.0 242956.0 146776.0 145299.0
##  [120] 256461.0 157273.0 170293.0 251862.0 287815.0 173402.0 146769.0
##  [127] 161397.0 135773.0 151240.0 241824.0 178743.0 160909.0 140957.0
##  [134] 156545.0 143199.0 151745.0 229309.0 254861.0 159545.0 143645.0
##  [141] 145491.0 130661.0 179806.0 191479.0 213775.0 146674.0 136913.0
##  [148] 146058.0 137231.0 153106.0 270295.0 315054.0 202539.0 172155.0
##  [155] 175498.0 147862.0 163339.0 221198.0 260549.0 161849.0 126401.0
##  [162] 154025.0 129388.0 153357.0 243271.0 211324.0 178725.0 150466.0
##  [169] 155214.0 136077.0 147621.0 209907.0 250089.0 156854.0 138015.0
##  [176] 142935.0 119825.0 148572.0 245232.0 268052.0 205679.0 181719.0
##  [183] 181846.0 146240.0 162050.0 224653.0 241187.0 167822.0 141820.0
##  [190] 146395.0 126275.0 141722.0 201631.0 211743.0 164808.0 158885.0
##  [197] 154253.0 128791.0 143851.0 205156.0 210765.0 157474.0 145577.0
##  [204] 144869.0 133297.0 143310.0 196669.0 207425.0 151433.0 145306.0
##  [211] 169011.0 186794.0 182056.0 244496.0 250045.0 187292.0 163162.0
##  [218] 169129.0 135575.0 154393.0 196019.0 209160.0 157230.0 144993.0
##  [225] 160384.0 143687.0 167387.0 211982.0 224815.0 164780.0 156473.0
##  [232] 161415.0 140135.0 153889.0 194942.0 214817.0 159968.0 143405.0
##  [239] 153890.0 135897.0 164717.0 254021.0 319865.0 217505.0 185584.0
##  [246] 189817.0 151754.0 152874.0 249262.0 283542.0 179703.0 139585.0
##  [253] 156259.0 131519.0 149029.0 229458.0 265735.0 178115.0 147840.0
##  [260] 151382.0 127872.0 153822.0 224371.0 247655.0 167229.0 140217.0
##  [267] 149795.0 124306.0 147889.0 231180.0 263403.0 170657.0 208268.0
##  [274] 201452.0 161280.0 175207.0 261633.0 284783.0 176256.0 146918.0
##  [281] 152104.0 134189.0 177986.0 196562.0 231014.0 158182.0 138236.0
##  [288] 171860.0 150942.0 147755.0 237509.0 266255.0 161032.0 138133.0
##  [295] 155821.0 130150.0 146654.0 215527.0 253798.0 160247.0 136578.0
##  [302] 152108.0 138288.0 188593.0 273130.0 307064.0 197117.0 176037.0
##  [309] 183333.0 162891.0 159996.0 244540.0 292804.0 170425.0 187991.0
##  [316] 159182.0 129556.0 161369.0 250917.0 290037.0 177915.0 152555.0
##  [323] 168996.0 135370.0 152956.0 229608.0 273432.0 169452.0 147986.0
##  [330] 167046.0 136741.0 160183.0 257947.0 354425.0 248377.0 203496.0
##  [337] 206228.0 186719.0 209627.0 251960.0 286953.0 206963.0 196903.0
##  [344] 208551.0 183048.0 199962.0 264214.0 319643.0 247739.0 224615.0
##  [351] 242840.0 221721.0 245301.0 299059.0 324394.0 343779.0 265132.0
##  [358] 201630.0 198375.0 227063.0 217157.0 273044.0 208788.0 405465.9
##  [365] 347068.5 496668.3 521388.7 315995.8 283933.6 351595.1 254783.5
##  [372] 282245.9 406345.6 456475.2 285943.1 249641.9 310704.2 252931.6
##  [379] 277379.3 411152.5 465501.8 296496.2 269739.5 303422.6 234063.5
##  [386] 264844.9 380838.3 439695.7 272729.9 236802.7 272548.7 225569.0
##  [393] 282855.7 303863.0 299998.0 214985.0 180542.0 190270.0 162112.0
##  [400] 175717.0 240940.0 249139.0 182044.0 162010.0 182976.0 157947.0
##  [407] 160134.0 255147.0 299729.0 193050.0 168587.0 171481.0 145280.0
##  [414] 166807.0 255773.0 235056.0 186126.0 159231.0 162456.0 144613.0
##  [421] 184072.0 429480.6 334158.9 325310.9 430482.7 378831.2 298411.0
##  [428] 294728.2 385211.9 501261.0 297661.1 271261.3 310927.8 233981.8
##  [435] 259745.2 370874.4 525402.9 313746.1 266588.0 319741.1 250400.0
##  [442] 276020.4 394767.1 458738.1 294048.3 273352.1 307726.6 232032.9
##  [449] 269390.9 426756.2 480275.5 329325.6 227833.0 222145.0 175226.0
##  [456] 190775.0 297111.9 319278.5 200069.7 172838.9 180832.5 152052.8
##  [463] 171982.6 253300.0 270253.0 176795.9 165316.0 197775.3 178496.5
##  [470] 186157.1 227965.4 262305.4 184750.8 171423.1 173458.9 143650.3
##  [477] 157121.3 233352.2 281722.1 166706.8 145001.0 163907.6 252978.2
##  [484] 198823.3 285064.4 315292.0 196322.7 174332.8 181574.9 148487.8
##  [491] 169429.8 268096.0 210572.8 178959.8 161410.6 165654.5 152593.9
##  [498] 173778.8 253342.1 310383.6 180361.9 158563.1 166948.0 141012.7
##  [505] 155391.7 237012.8 275004.8 172440.1 147170.0 156817.3 137604.8
##  [512] 167732.4 270734.1 352715.2 232747.4 199679.8 195493.2 158628.3
##  [519] 173425.6 255098.2 298182.8 184713.1 153929.6 171354.5 139568.4
##  [526] 175078.2 288595.0 196388.7 203934.0 171726.6 174330.3 148041.5
##  [533] 158803.4 247238.3 282880.9 193543.4 151639.3 159011.2 140099.9
##  [540] 163984.3 233325.2 270481.1 193755.8 347219.6 390813.8 288799.8
##  [547] 289130.7 403329.5 477001.9 317146.1 264221.8 301686.7 238018.4
##  [554] 267158.8 379943.7 369119.2 294631.2 271954.6 332075.3 245724.5
##  [561] 266066.2 364000.9 396616.6 297509.6 257608.2 304069.1 237145.7
##  [568] 287240.5 342228.9 383901.8 291216.4 256697.0 307711.2 258139.3
##  [575] 227181.3 285339.2 293770.8 216069.1 189204.5 185716.1 158289.6
##  [582] 175001.3 225347.2 238107.3 183748.5 154682.6 175462.6 154187.7
##  [589] 181224.0 256980.3 280854.6 199989.6 172931.4 185372.8 161694.6
##  [596] 172877.6 230664.8 251529.2 191927.7 159008.3 174562.1 148645.6
##  [603] 184555.5 279284.6 328072.9 419468.4 361543.3 399488.0 274916.5
##  [610] 306029.6 456119.4 542624.6 332527.0 281738.7 334461.0 236246.6
##  [617] 285118.8 412331.3 483844.2 322557.8 306315.5 339881.3 234539.7
##  [624] 270423.6 424190.3 489017.9 301604.6 281400.2 325418.7 232100.3
##  [631] 255826.3 415323.6 488492.6 300600.1 297270.7 421319.3 314106.9
##  [638] 337466.0 488112.0 542030.3 319723.5 290411.3 337023.6 242559.6
##  [645] 350580.5 381063.6 439153.9 308077.4 273305.9 338017.4 267781.0
##  [652] 282585.1 426339.1 502293.3 320542.4 268577.0 317293.7 241560.4
##  [659] 261200.6 392114.7 476067.3 284784.4 249220.1 300899.3 222880.2
##  [666] 275187.4 424910.9 363067.7 519232.5 343948.0 382202.5 279404.3
##  [673] 291440.7 424269.5 498982.3 298188.1 276371.8 315737.8 229578.5
##  [680] 275051.6 405434.1 514483.2 320294.7 278042.8 307401.7 253722.4
##  [687] 285957.5 391566.1 482586.1 315183.5 273281.1 313890.9 235714.3
##  [694] 281230.0 432606.1 543247.1 428088.8 368561.1 382582.5 280652.1
##  [701] 320977.5 417563.1 546761.4 391494.9 348685.8 410567.5 291936.9
##  [708] 317018.1 434492.2 522822.8 425669.7 420466.8 406694.0 360477.1
##  [715] 364679.0 449397.1 529028.3 490758.2 475355.7 449102.3 401854.2
##  [722] 337823.8 346707.5 359718.0 413535.7 370377.7 282875.0 284799.7
##  [729] 352874.4 267312.8 215502.1 230227.6 182385.9 190076.9 261700.9
##  [736] 331613.0 197527.7 172916.2 181581.3 160873.2 183475.5 288426.3
##  [743] 325833.1 203409.9 178523.5 189735.0 157668.7 178822.2 269758.8
##  [750] 296561.9 180873.0 155998.0 171530.1 146541.0 189210.8 301041.6
##  [757] 389954.2 239899.4 208233.1 207420.6 176824.0 190362.6 272478.0
##  [764] 287296.1 198896.2 169499.6 172341.5 159941.7 207827.4 226863.5
##  [771] 189035.0 192378.6 236416.4 214186.6 181948.6 212016.8 262266.9
##  [778] 338207.0 196747.3 187260.0 179290.5 152255.6 189490.8 296928.6
##  [785] 402789.0 248840.1 221231.7 221144.2 186864.9 200771.4 303724.0
##  [792] 329978.0 221274.3 209066.2 196797.6 167976.6 189363.3 287007.8
##  [799] 330108.3 211844.0 184646.6 181460.4 164252.0 210532.8 297590.2
##  [806] 320583.5 211516.4 180545.6 189016.3 159054.8 194418.7 273913.7
##  [813] 317237.2 200597.7 199211.7 260693.9 231146.4 240140.4 299513.5
##  [820] 324691.2 238320.7 211640.9 214879.6 173543.1 195900.8 285422.9
##  [827] 328175.6 201698.2 182022.2 206746.8 189224.0 203730.5 299710.6
##  [834] 348888.9 228836.2 195669.8 197303.5 177838.0 193409.8 275928.0
##  [841] 312118.9 195023.3 170086.6 184228.0 190483.5 309666.1 325500.4
##  [848] 362793.3 254782.0 227581.4 229586.9 190300.9 212320.9 326621.7
##  [855] 256707.6 239950.6 204035.0 208834.3 172561.8 222143.8 338783.0
##  [862] 392400.0 231566.6 205565.5 219039.1 188560.9 204723.8 305234.6
##  [869] 360877.4 253845.0 218470.2 234518.3 196568.0 256271.5 407567.2
##  [876] 511203.3 417122.8 375666.8 403911.0 301352.7 328171.9 464018.7
##  [883] 569749.0 354295.7 313259.8 341191.3 248065.7 303813.5 430946.2
##  [890] 513090.9 338502.2 342322.4 364956.1 272748.7 335740.1 490718.8
##  [897] 405074.9 344391.2 323181.9 350053.0 258099.6 298247.6 430996.1
##  [904] 528912.4 324198.2 326055.6 470543.3 356507.7 399224.4 533329.9
##  [911] 519888.9 359720.1 306808.9 399941.7 312194.1 322732.4 405364.1
##  [918] 469498.2 323429.4 302241.4 355368.0 293701.0 325329.8 422336.3
##  [925] 468922.9 339664.7 299832.6 333165.1 263786.6 301333.5 381750.4
##  [932] 440293.9 317205.5 282849.4 324089.6 263595.8 328760.2 511251.1
##  [939] 524731.8 391042.7 349915.9 368385.0 276042.3 311399.6 370916.7
##  [946] 323678.3 392416.0 329890.1 349608.8 277403.0 385882.1 492820.1
##  [953] 473570.8 348131.4 307598.5 338124.7 276466.0 317804.5 393233.1
##  [960] 437910.6 334308.1 310822.3 336598.9 268717.6 316715.5 444344.5
##  [967] 513127.5 389380.4 451319.1 464525.4 358435.9 372872.2 545626.9
##  [974] 616773.6 384504.7 328000.0 369268.1 280890.7 322584.5 472754.5
##  [981] 546711.5 337937.5 334826.6 395899.4 290574.6 319025.2 480835.8
##  [988] 559209.2 357203.0 314663.2 360131.1 274412.9 319862.6 461417.5
##  [995] 545635.3 332627.5 298849.7 375304.8 395489.9 412480.0 605299.6
## [1002] 667910.5 409239.9 340063.9 375993.1 281981.0 424475.2 423698.8
## [1009] 495042.7 347861.5 267742.3 357180.5 281102.5 336445.2 493557.7
## [1016] 559176.8 348644.1 318378.7 354599.6 268801.3 317892.0 442228.0
## [1023] 529488.4 317842.4 293573.3 314592.6 244752.9 303908.1 441528.1
## [1030] 397013.6 386927.5 547695.9 440339.8 323681.6 329357.5 483088.5
## [1037] 579957.8 336585.7 309336.5 387229.8 244509.5 307238.5 449429.2
## [1044] 566775.4 362567.3 311767.8 363735.4 282805.6 332294.0 450939.1
## [1051] 545786.9 360896.0 305097.7 358997.7 267631.9 316686.0 437186.2
## [1058] 580988.4 395547.6 432504.9 465334.0 333036.3 377201.3 485246.5
## [1065] 614647.1 433837.8 359011.9 404125.9 328620.5 351776.6 478540.5
## [1072] 568795.8 387337.7 378103.6 436770.3 384431.4 435312.4 528456.7
## [1079] 581448.5 499726.2 520204.0 568770.5 497713.0 437337.6 413311.9
## [1086] 365548.8 398811.7 490101.7 427365.7 526644.0 660848.1 494334.2
## [1093] 421425.3 440002.8 324676.7 353317.8 514767.8 585899.8 348082.5
## [1100] 316359.5 347298.4 277801.4 327304.0 494916.6 563633.9 358244.1
## [1107] 309266.0 395699.2 272574.7 320448.5 465111.2 523837.3 332162.8
## [1114] 301251.9 342945.4 263094.8 319881.7 503629.1 587306.6 434383.9
## [1121] 366993.5 394391.6 298919.4 353269.8 437250.1 330186.5 324413.2
## [1128] 409769.3 389318.9 307548.0 338633.0 437553.7 504309.0 382447.2
## [1135] 364012.8 388859.0 276005.5 333130.8 457357.9 540335.2 350412.3
## [1142] 320433.1 355730.8 272377.1 314854.4 471167.4 538887.1 368794.2
## [1149] 440298.1 450724.5 329458.0 363172.3 534252.2 588029.7 367665.6
## [1156] 316725.6 364907.2 279582.0 326087.9 462500.3 518177.1 331437.4
## [1163] 328944.2 393335.0 288876.7 331396.5 455968.1 500560.6 343100.0
## [1170] 295950.7 368277.2 314574.1 313116.7 397462.5 447119.3 332162.5
## [1177] 306377.0 348440.5 311121.3 417783.5 594849.7 692879.2 400124.5
## [1184] 354813.1 390764.1 291944.7 339236.6 478171.6 553962.9 337192.6
## [1191] 328033.4 356020.4 268332.1 325316.2 439220.3 741244.1 774111.4
## [1198] 632660.9 585920.0 597476.7 454771.3 538819.3 557393.3 338659.0
## [1205] 377142.6 356345.7 269901.9 349303.0 545553.9 665095.8 542026.2
## [1212] 387182.5 414267.2 312126.5 373266.3 544143.6 406535.5 381160.0
## [1219] 328716.0 383262.6 307095.3 325211.3 473743.3 557182.2 371337.1
## [1226] 321526.8 365997.2 284173.0 346338.6 461155.5 534171.0 348836.7
## [1233] 289793.5 331838.4 281615.8 423390.9 397858.3 495796.6 383925.8
## [1240] 486524.9 416423.1 340132.9 360083.9 505333.4 582898.0 365706.5
## [1247] 321906.4 362524.5 269844.8 314176.2 425214.5 486396.1 353006.1
## [1254] 303267.0 372069.8 304986.7 383563.8 497934.6 420956.3 369478.7
## [1261] 333177.9 353238.4 273641.4 316357.7 438390.8 492091.4 355404.8
## [1268] 292903.1 330505.6 269992.9 411626.7 525960.0 576086.7 433759.1
## [1275] 348410.8 388027.0 286932.1 336089.9 431364.9 471440.6 337019.4
## [1282] 302502.2 334369.1 261372.1 326367.7 439766.2 474843.3 343807.2
## [1289] 306656.2 354317.0 288828.2 313986.1 399392.0 425174.4 339340.5
## [1296] 294460.8 339372.5 260636.4 322076.2 446431.5 504993.6 434916.1
## [1303] 370435.2 399429.0 302264.8 339102.6 421584.7 458621.2 349854.2
## [1310] 300570.5 337207.5 284272.5 358675.7 340230.1 382917.3 348033.6
## [1317] 352786.3 373172.0 277207.4 329083.0 417951.1 458306.3 354765.6
## [1324] 319851.0 346442.2 274231.0 315997.8 391514.7 443223.5 331724.6
## [1331] 304524.6 381155.8 338440.7 418751.1 536551.0 574919.7 408240.6
## [1338] 365935.5 393172.8 297387.6 327585.3 470524.8 544330.7 335844.4
## [1345] 291627.5 333137.2 277469.8 340069.5 466964.5 539694.5 343172.2
## [1352] 322120.7 353503.0 262546.8 326668.3 445565.9 508538.4 331042.2
## [1359] 288279.6 334374.9 256807.9 347445.3 608232.4 685731.8 422783.9
## [1366] 356339.1 401434.6 284469.5 351032.3 483704.2 583564.5 358195.2
## [1373] 303556.3 350186.3 269374.8 331432.5 452758.0 593811.0 363633.8
## [1380] 309542.5 341612.0 292229.4 339759.5 438222.5 549176.6 330197.3
## [1387] 303713.4 336867.5 265582.4 316608.9 426465.6 506636.8 318896.6
## [1394] 386140.3 481823.0 380573.0 418430.9 446728.2 607668.7 419944.2
## [1401] 384282.4 415179.7 293912.1 347040.6 481531.7 639242.1 372847.4
## [1408] 338379.2 441808.9 315953.4 373908.4 481869.2 608879.4 385165.8
## [1415] 330607.5 395565.0 292831.4 358721.3 465100.0 566030.3 390483.5
## [1422] 318030.1 413701.1 423594.3 404948.2 557976.7 674883.4 444185.2
## [1429] 510508.0 462025.6 345984.6 392411.3 475630.0 595882.7 410557.5
## [1436] 368920.4 443581.8 369755.4 424198.1 537535.9 647831.7 501520.8
## [1443] 522373.7 573738.0 525120.5 601798.1 592323.4 513937.1 392031.5
## [1450] 452519.9 411588.8 544195.0 496313.4 737238.0 559340.8 515173.8
## [1457] 395948.5 430914.7 574376.1 639935.5 402034.7 363447.5 408179.4
## [1464] 312713.6 384748.9 532191.5 640421.0 418809.1 381744.3 414257.3
## [1471] 300270.3 359032.5 564619.4 627341.6 404989.5 356288.5 385545.2
## [1478] 282025.5 357593.3 494317.6 563023.1 344587.2 385595.7 497565.9
## [1485] 368365.8 410262.2 577135.9 600102.3 394942.3 366208.8 410151.3
## [1492] 312077.0 377972.8 479058.9 573369.8 400505.9 312221.4 384814.9
## [1499] 323944.9 379077.2 561712.3 438901.4 373573.1 354786.7 375507.2
## [1506] 322535.4 376564.3 476880.6 376740.9 381927.2 520136.1 499310.2
## [1513] 407898.1 425594.8 576151.6 645301.2 390805.6 380877.3 404202.2
## [1520] 320176.0 388727.0 516446.9 627965.7 384266.3 346121.6 429516.8
## [1527] 348306.5 406969.3 540616.5 638106.7 419307.9 375820.1 404179.5
## [1534] 341858.8 390287.1 535592.1 594608.3 395793.4 340583.0 428750.4
## [1541] 314090.4 444690.8 782010.7 601792.8 397390.5 411686.8 442585.3
## [1548] 340652.0 372748.8 512159.5 597139.7 365009.3 376358.7 407646.7
## [1555] 367009.4 360725.0 447295.7 529100.7 386402.2 353225.3 388741.9
## [1562] 323269.3 381320.3 483353.6 593372.5 358139.5 344886.7 386422.4
## [1569] 296278.1 372999.6 512914.9 508631.8 672300.3 439158.9 473210.9
## [1576] 353384.8 378791.7 516483.4 625263.6 385273.4 346126.8 390585.0
## [1583] 307280.2 363369.5 559962.3 431840.7 400089.2 393390.9 425287.0
## [1590] 304206.8 369917.3 521946.9 613917.6 375577.6 347380.3 375145.9
## [1597] 289416.1 464841.2 439138.2 534458.1 369374.3 355445.3 420619.5
## [1604] 375880.0 445102.8 590609.0 716517.2 441056.5 384463.4 428390.5
## [1611] 315565.5 378511.5 504078.4 630622.2 368640.6 354659.4 391660.2
## [1618] 320264.5 408481.4 551790.6 471786.8 389198.4 387088.0 390078.4
## [1625] 295067.5 363312.5 483357.3 557800.6 356927.8 325677.7 368282.0
## [1632] 296359.6 386850.5 590306.1 667458.7 453488.7 397778.2 424935.0
## [1639] 333543.9 377091.0 492070.9 549406.1 386486.3 347868.1 369429.6
## [1646] 297521.3 353842.2 469016.4 552144.7 389266.2 343912.1 385352.3
## [1653] 327844.1 362561.2 451583.5 500142.6 385172.0 331758.9 361277.0
## 
## Labels: 
## [1] "Level 0" "Level 1" "Level 2" "Level 3"

val.hts = window(train.hts, start = 2017.542) # 20 dates
summary(val.hts)

## Hierarchical Time Series 
## 4 Levels 
## Number of nodes at each level: 1 2 9 45 
## Total number of series: 57 
## Number of observations per series: 20 
## Top level series: 
## Time Series:
## Start = 2017.54209445585 
## End = 2017.59411362081 
## Frequency = 365.25 
##  [1] 307073.0 383286.5 492941.0 556749.5 421266.4 462371.9 464175.3
##  [8] 341257.8 376220.7 467873.3 520594.2 370774.0 328836.5 352419.5
## [15] 292493.3 372508.4 368959.1 424217.0 353314.0 356006.0
## 
## Labels: 
## [1] "Level 0" "Level 1" "Level 2" "Level 3"

SMAPE Function

smape = function(y_pred, y_actual) {
  error = mean(abs(y_pred - y_actual)/(abs(y_actual) + abs(y_pred)))
  return(error)
}

Base Model

rw.top.down = forecast(
  data.hts, h = 20, method = "tdfp", fmethod = "rw",
  keep.fitted = TRUE, keep.resid = TRUE
)

rw.bottom.up = forecast(
  data.hts, h = 20, method = "bu", fmethod = "rw",
  keep.fitted = TRUE, keep.resid = TRUE
)

rw.middle.out = forecast(
  data.hts, h = 20, method = "mo", fmethod = "rw", level = 2,
  keep.fitted = TRUE, keep.resid = TRUE
)

rw.comb = forecast(
  data.hts, h = 20, method = "comb", fmethod = "rw",
  keep.fitted = TRUE, keep.resid = TRUE
)

base_results = rbind(smape(rw.top.down$bts, val.hts$bts),
                smape(rw.bottom.up$bts, val.hts$bts),
                smape(rw.middle.out$bts, val.hts$bts),
                smape(rw.comb$bts, val.hts$bts))
rownames(base_results) = c("top.down",  "bottom.up",  "middle.out", "combination")
base_results

##                  [,1]
## top.down    0.1521026
## bottom.up   0.1521026
## middle.out  0.1521026
## combination 0.1521026

Irrespective of method, hierarchical time series with random walk gives the same SMAPE of 0.152.

Base model plots (Top Down)

rw.top.down$bts[,1:5] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

val.hts$bts[,1:5] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

rw.top.down$bts[,26:30] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,87000)

val.hts$bts[,26:30] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,87000)

Forecast using different methods

top.down = forecast(
  data.hts, h = 20, method = "tdfp", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE
)

bottom.up = forecast(
  data.hts, h = 20, method = "bu", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE
)

middle.out = forecast(
  data.hts, h = 20, method = "mo", fmethod = "arima", level = 2,
  keep.fitted = TRUE, keep.resid = TRUE
)

comb = forecast(
  data.hts, h = 20, method = "comb", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE
)

# More combination variations were tried
# comb_wls = forecast(
#   data.hts, h = 20, method = "comb", fmethod = "arima", weights = "wls",
#   keep.fitted = TRUE, keep.resid = TRUE
# )
#  
# comb_wls_alg = forecast(
#   data.hts, h = 20, method = "comb", fmethod = "arima", weights = "wls", algorithms = "lu",
#   keep.fitted = TRUE, keep.resid = TRUE
# )
# 
# comb_nseries = forecast(
#   data.hts, h = 20, method = "comb", fmethod = "arima", weights = "nseries",
#   keep.fitted = TRUE, keep.resid = TRUE
# )

Compare SMAPE Errors

results = rbind(smape(top.down$bts, val.hts$bts),
                smape(bottom.up$bts, val.hts$bts),
                smape(middle.out$bts, val.hts$bts),
                smape(comb$bts, val.hts$bts))
rownames(results) = c("top.down",  "bottom.up",  "middle.out", "combination")
results

##                   [,1]
## top.down    0.10123404
## bottom.up   0.09917318
## middle.out  0.09957984
## combination 0.09841618

Forecast Plots - Arima

comb$bts[,1:5] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

val.hts$bts[,1:5] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

comb$bts[,26:30] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,90000)

val.hts$bts[,26:30] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,90000)

Forecast - Arima with XREG

Oil Prices

oil_price <- unique(train[,c("date","price")])
train_oil_price <- oil_price$price[1:1659]
test_oil_price <- oil_price$price[1660:1679]

comb_oil_xreg = forecast(
  data.hts, h = 20, method = "comb", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE, xreg = train_oil_price, newxreg = test_oil_price
)

smape(comb_oil_xreg$bts, val.hts$bts)

## [1] 0.09532138

Oil Price, National Flag, Holiday Flag and Pay Day Flag

x_var <- unique(train[,c(2,7,10:12)])
train_x_var <- as.matrix(x_var[1:1659,2:5])
test_x_var <- as.matrix(x_var[1660:1679,2:5])

comb_all_xreg = forecast(
  data.hts, h = 20, method = "comb", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE, xreg = train_x_var, newxreg = test_x_var
)

smape(comb_all_xreg$bts, val.hts$bts)

## [1] 0.09712447

Forecast Plots - Arima +Xreg

comb_all_xreg$bts[,1:5] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

val.hts$bts[,1:5] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

comb_all_xreg$bts[,26:30] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,90000)

val.hts$bts[,26:30] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Quinto Type A")+
  ylim(8000,90000)

Oil Price, National Flag, Holiday Flag, Pay Day Flag, Perishable and Promotion

# total n_promotions in unique day
total_n_promotion = train %>% dplyr::group_by(date) %>% dplyr::summarise(total_n_promotion =sum(n_promotion))
## total number of perishable
total_perishable = train %>% dplyr::group_by(date) %>% dplyr::summarise(total_perishable = sum(perishable))

x_var$total_n_promotion <- total_n_promotion$total_n_promotion
x_var$total_perishable <- total_perishable$total_perishable

train_x_var_all <- as.matrix(x_var[1:1659,2:7])
test_x_var_all <- as.matrix(x_var[1660:1679,2:7])

comb_all_xreg_n = forecast(
  data.hts, h = 20, method = "comb", fmethod = "arima",
  keep.fitted = TRUE, keep.resid = TRUE, xreg = train_x_var_all, newxreg = test_x_var_all
)

# Overall SMAPE
smape(comb_all_xreg_n$bts, val.hts$bts)

## [1] 0.08731203

SMAPE Across each group

smape.values <- data.frame()
for (i in 1:45){
  smape.values<-rbind(smape.values, smape(comb_all_xreg_n$bts[,i], val.hts$bts[,i]))
}

smape.values

##    X0.0688789317414289
## 1           0.06887893
## 2           0.06268975
## 3           0.06502918
## 4           0.08112434
## 5           0.06468529
## 6           0.07398567
## 7           0.06344047
## 8           0.10173905
## 9           0.08843921
## 10          0.05633753
## 11          0.12220925
## 12          0.21415143
## 13          0.07650701
## 14          0.10987188
## 15          0.10792249
## 16          0.06039761
## 17          0.06796901
## 18          0.04701139
## 19          0.08949439
## 20          0.07445672
## 21          0.09936316
## 22          0.14272833
## 23          0.11989348
## 24          0.11642535
## 25          0.13296809
## 26          0.08644438
## 27          0.09366623
## 28          0.04197678
## 29          0.09213388
## 30          0.05533820
## 31          0.06095472
## 32          0.17208388
## 33          0.06409111
## 34          0.07850777
## 35          0.14054259
## 36          0.04562756
## 37          0.13680483
## 38          0.07543752
## 39          0.06395193
## 40          0.06578375
## 41          0.09326812
## 42          0.08644490
## 43          0.03162136
## 44          0.05777897
## 45          0.07886371

smape.final <- colSums(smape.values)/45; smape.final

## X0.0688789317414289 
##          0.08731203

Forecast Plots - Arima Xreg

comb_all_xreg_n$bts[,1:5] %>%
  autoplot(facet=TRUE) +
  xlab("Year") + ylab("Total Sales") + ggtitle("Forecasted Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

val.hts$bts[,1:5] %>%
  autoplot(facet=TRUE, color='red')+
  xlab("Year") + ylab("Total Sales") + ggtitle("Actual Grocery Sales by City, Type & Item Family - Guayaquil Type A")+
  ylim(800,10000)

Cross-Validation Function using Arima & Xreg

hts.CV = function(y, window, h) {
  cv.error = matrix(ncol = 1)
  
  for (i in seq(1, nrow(y) - (window + h) + 1)) {
    train_start = i
    train_end = i + window - 1
    test_start = i + window
    test_end = i + window + h - 1
   
    train_y = subset(y, start = train_start, end = train_end)
    test_y = subset(y, start = test_start, end = test_end)

    train_hts = hts(train_y, characters = c(2, 1, 5))
    test_hts = hts(test_y, characters = c(2, 1, 5))
    
    train_x_var_all_cv <- as.matrix(x_var[train_start:train_end,2:7])
    test_x_var_all_cv <- as.matrix(x_var[test_start:test_end,2:7])
    
    fc = forecast(train_hts, h = h, method = "comb", fmethod = "arima",
    keep.fitted = TRUE, keep.resid = TRUE, xreg = train_x_var_all_cv, newxreg = test_x_var_all_cv)
    print(paste("Done", as.character(i)))
    acc = smape(fc$bts, test_hts$bts)
    cv.error = rbind(cv.error, acc)
  }
  
  return(cv.error)
}

CV Using 10 Folds

start.time = Sys.time()
cv.acc = hts.CV(y = train.ts, window = 1650, h = 20) # takes 11 mins to run one fold

## [1] "Done 1"
## [1] "Done 2"
## [1] "Done 3"
## [1] "Done 4"
## [1] "Done 5"
## [1] "Done 6"
## [1] "Done 7"
## [1] "Done 8"
## [1] "Done 9"
## [1] "Done 10"

end.time = Sys.time()
(time.taken = end.time - start.time)

## Time difference of 6.3386 hours

Due to limited computational power, only 3-fold CV was implemented.

CV Errors

# Mean SMAPE
mean(cv.acc, na.rm = T)

## [1] 0.08551528

# Standard Deviation of SMAPE
sd(cv.acc, na.rm = T)

## [1] 0.004075772